n_hub = 500
k = 100
hub_scores <- hubness$score[hubness$Dimension==(cell_nb-1) & hubness$p==2 & hubness$k==k]
names(hub_scores) <- colnames(seurat)
k_neighbors <- colnames(seurat)[unique(unlist(as.list(neighbor_list[order(hub_scores, decreasing  =T)[1:n_hub],])))]
seurat$connection <- ifelse(colnames(seurat) %in% k_neighbors, "connected", "free")

print(paste0(length(k_neighbors),
             " cells are connected to hubs, i.e. ",
             floor(length(k_neighbors)/ncol(seurat)*100),
             "% of the total cells"))
## [1] "5762 cells are connected to hubs, i.e. 63% of the total cells"
DimPlot(seurat, group.by = "connection", reduction = "umap")

# nFeature/nCount
ggarrange(FeaturePlot(seurat,"nFeature_RNA"),
          FeaturePlot(seurat,"nCount_RNA"))

ggarrange(vlnplot_hub_spec_connect(seurat,"nFeature_RNA","connection",800,8000,7900),
          vlnplot_hub_spec_connect(seurat, "nCount_RNA", "connection",0,3e6,3e6))
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

ggarrange(histoplot_hub_spec(seurat, seurat$nFeature_RNA, "nFeature", seurat$connection),
          histoplot_hub_spec(seurat, seurat$nCount_RNA, "nCount", seurat$connection))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# Dropout rate
FeaturePlot(seurat,"dropout")

vlnplot_hub_spec_connect(seurat,"dropout","connection",0.8,1.05,1.02)
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Warning: Removed 485 rows containing non-finite values (stat_ydensity).
## Warning: Removed 485 rows containing non-finite values (stat_signif).
## Warning: Removed 485 rows containing missing values (geom_point).

histoplot_hub_spec(seurat, seurat$dropout, "dropout", seurat$connection)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# Percentage of mito/ribo
ggarrange(FeaturePlot(seurat,"percent.mt"),
          FeaturePlot(seurat,"percent.rb"))
## Warning in FeaturePlot(seurat, "percent.mt"): All cells have the same value (0)
## of percent.mt.

ggarrange(vlnplot_hub_spec_connect(seurat, "percent.mt", "connection",0,5,4),
          vlnplot_hub_spec_connect(seurat, "percent.rb", "connection",0,21,20))
## Warning in SingleExIPlot(type = type, data = data[, x, drop = FALSE], idents =
## idents, : All cells have the same value of percent.mt.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

ggarrange(histoplot_hub_spec(seurat, seurat$percent.mt, "percent.mt", seurat$connection),
          histoplot_hub_spec(seurat, seurat$percent.rb, "percent.rb", seurat$connection))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# scEntropy 
FeaturePlot(seurat, reduction = "umap", features = "scentropy")

vlnplot_hub_spec(seurat, "scentropy", "connection",6,6.9,6.8)
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Removed 1 rows containing missing values (geom_hline).

histoplot_hub_spec(seurat, seurat$scentropy, "scentropy", seurat$connection)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# stemID
FeaturePlot(seurat, reduction = "umap", features = "stemID")

vlnplot_hub_spec(seurat, "stemID", "connection",0.5,0.8,0.78)
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Removed 1 rows containing missing values (geom_hline).

histoplot_hub_spec(seurat, seurat$stemID, "stemID", seurat$connection)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.